med_pvals <- function(sim.name){
  # Initialize an empty list to store p-values from all 10 simulations
  all_p_values <- list()
  
  # Loop through each simulation result and extract the p-values
  for (i in 1:10) {
    p_values <- c(
      sim.name[[i]]$d.avg.p,
      sim.name[[i]]$z.avg.p,
      sim.name[[i]]$n.avg.p
    )
    
    # Store the extracted p-values in the list
    all_p_values[[i]] <- p_values
}
  # Convert list of p-values to a matrix (for easier manipulation)
  all_p_values_matrix <- do.call(rbind, all_p_values)
  # Assign column names to the p-values (so you know which term each corresponds to)
  colnames(all_p_values_matrix) <- c(
    "ACME change p-value", 
    "ADE change p-value", 
    "Total change p-value"
  )
  
  # Calculate the average p-value for each of the 10 values across all simulations
  avg_p_values <- colMeans(all_p_values_matrix)
  # Print the average p-values
  avg_p_values
}